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■^ ' Abstract 

^^ ■ We propose a two dimensional (2D) adaptive nodes technique for ir- 

regular regions. The method is based on equi-distribution principal and 
dimension reduction. The mesh generation is carried out by first pro- 
2 ' ducing some adaptive nodes in a rectangle based on equi-distribution 

along the coordinate axes and then transforming the generated nodes 
to the physical domain. Since the produced mesh is applied to the 
meshless-type methods, the connectivity of the points is not used and 
only the grid points are important, though the grid lines are utilized in 
the adapting process. The performance of the adaptive points is exam- 
\^^ ' ined by considering a collocation meshless method which is based on 

\^ , interpolation in terms of a set of radial basis functions. A generalized 

^S| ■ thin plate spline with sufficient smoothness is used as a basis function 

and the arc-length is employed as a monitor in the equi-distribution 
process. Some experimental results will be presented to illustrate the 
effectiveness of the proposed method. 

Key words. Adaptive mesh, Equi-distribution, Irregular regions, Colloca- 
tion meshless method, Dimension reduction method. 



'j_j ■ 1 Introduction 



Mesh free methods are now well known in the numerical solution of partial 
differential equations (PDEs). The most attractive feature of these tech- 
niques is that discretization of the domain or boundary is not required. This 
feature considerably reduces the computational complexity of the method. 
While the meshless techniques are often divided into two categories: bound- 
ary (see for example [20]) and domain type methods, the current work con- 
cerns with the latter category. 
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The use of radial basis functions (RBFs) for solving PDEs, first presented by 
Kansa, is a fully mesh free approach and falls into the domain type matheds 
[121124]. This method can be easily applied to the case of higher dimensional 
spaces due to the nature of the RBFs. Despite a good performance of RBFs 
in approximating multi-variate functions, they involve ill-conditioning, es- 
pecially for large scale problems. Another difficulty concerns their compu- 
tational efficiency, due to the dense matrices arising from interpolation. To 
tackle the above difhculties some sort of localization, such as domain decom- 
position methods (DDM) [8] and compactly supported RBFs (CS-RBFs) [5], 
the most important of which was introduced by Wendland [21], have been 
recommended. 

In the DDM the domain is divided into some subdomains and the PDE is 
solved for each subproblem followed by assembling the global solution. As 
a result, the ill-conditioning is avoided and the computational efficiency is 
improved due to working with small size matrices. On the other hand, using 
the CS-RBFs results in sparse matrices, which again improve the condition- 
ing and computational efficiency of the method. 

This work involves a different approach which can still be used with the 
above proposed methods. As was highlighted before, in using the classi- 
cal RBFs, increasing the size of the problem itself affects the conditioning. 
Consequently, reducing the number of nodes can improve the conditioning. 
One way to achieve this goal is to apply a set of adaptive nodes rather than 
uniform ones. As is well known, the main idea in adaptive meshes is to 
use a minimum number of nodes while still having the desired accuracy. 
This is achieved by allocating more mesh points to the areas where they are 
required. The adaptive mesh strategies often fall into two categories: the 
equi-distribution principle [7] and the variational principle [23j. The most 
popular technique, which has been widely used in the literature, is based on 
the equi-distribution strategy, which is also employed in this work. In this 
approach the mesh distribution is carried out in such a way that some mea- 
sure of error, called a monitor function, is equalized over each subinterval. 

Much effort has been devoted to generating adaptive meshes in two and 
three dimensional spaces, based on both equi-distribution and variational 
principles (see for example [18^ [2l [15]). In the literature two major meth- 
ods, namely transformation [6] and dimension reduction [19], have been 
employed to produce 2D meshes. The first category is based on mapping 
the physical domain into a simple domain with a uniform mesh and it leads 
to solving a differential equation in order to obtain an adaptive mesh. In 



the dimension reduction method, which is also employed in this paper, the 
equi-distribution process is reduced to a ID case. 

A method based on equi-distributing along the grid lines in the coordi- 
nate directions has been presented to produce mesh points in rectangular 
regions [16]. Also a generalization of this method to the case of three di- 
mensions was proposed in [17]. However, due to the use of grid lines in the 
coordinate directions, this method was limited to the case of rectangular and 
cubic domains, respectively, in the case of 2D and 3D. The purpose of the 
current work is to extend the above method to more general cases with irreg- 
ular boundaries in 2D. This is carried out by first, generating some adaptive 
nodes in a rectangle, then mapping the generated nodes to the physical do- 
main employing a suitable transformation. Of course, the mesh produced 
by the proposed method neither precisely satisfies the equi-distributing con- 
dition nor concerns about properties such as orthogonality, which are often 
required in numerical mesh-dependent methods. Instead, the mesh points 
are suitable for any meshless methods in which the mesh points, rather than 
mesh lines, are important. 

While a part of researches in the adaptive mesh community concerns 
with constructing monitor functions for different applications [3] , the current 
study focus on the mesh generation strategy using a well known monitor, 
namely arc-length j22j . In addition, we remark that the connectivity of the 
mesh are not used in this work, although they are used in the mesh points 
generation process. 

This paper is organized as follows. In section [2] the adaptive mesh technique 
in ID is reviewed. A generalization of this method to the case of 2D is 
discussed in section [3l The new mesh generation method is presented in 
section m In section [5] the collocation meshless method is reviewed. Some 
numerical results are given in section [6j 

2 Adaptive mesh 

We now introduce the concept of equi-distribution in the case of ID. 

Definition 1 (Equi-distributing) Let M he a non-negative piecewise con- 
tinuous function on [a, b] and c be a constant such that n = - f M{x)dx is 
an integer. The mesh 

11: a = xq < xi < ■ ■ ■ < Xn = b, 



is called equi- distributing (e.d.) on [a, 6] with respect to M and c if 
M{x)dx = c, J = 0, 1, • • • , n — 1. 

A suitable algorithm to produce an e.d. mesh has been given in |13j . 
In Definition [1] the function /, often called a monitor, is dependent on the 
solution of the underlying PDE and its derivatives. The arc-length monitor 

M = Vl + «^ (1) 

which is used in this work has been widely used in the literature (see for ex- 
ample [221 [T]). The function n in ([T| is the solution of the underlying PDE, 
X is the coordinate, in the direction of which the adaptivity is performed, 
and Ux is the partial derivative with respect to x. To find more details about 
the monitors for different applications see for example [U [3] . 



3 Adaptive nodes in a rectangle 

A natural extension of Definition [1] to the case of 2D is as follows. 

Definition 2 (2D Equi-distributing) Given a 2D domain VL, a 2D adap- 
tive mesh based on equi-distributing will be a mesh obtained by dividing the 
domain $7 into n subdomains Vli such that 

M{x,y)dxdy = constant, 

where M is a suitable monitor function. 

Obviously an infinite number of adaptive meshes based on Definition [2] exist. 
However, obtaining even one of these meshes can be a complicated process. 
A method based on dimension reduction was proposed for a rectangular re- 
gion in [16]. Since the current work is a development of the above mentioned 
method, here it is briefly reviewed. To do so, we start with a uniform mesh 
in a rectangle in the form 

{{x,y)\ ai<x<bi, a2<y < 62}- 

Let the underlying uniform mesh points be 

{{xij,yij)\ i = 0,l,--- ,iVi, j = 0,l--- ,iV2}, 



a 





c 



Figure 1: The three stages of the adaptive mesh method are shown in Fig- 
ures (a), (b) and (c) respectively and in each direction the grid curves are 
displayed with the quadrilateral formed. 



where 



Xi 



U 



and 



h. 



ai+i- hi, yij = yj = a2 + j ■ /i2, 

m = 1,2. 



N„ 



The equi-distribution process is performed in three stages. In the first stage, 
equi-distribution is performed in the direction of the x-axis. More precisely, 
for each horizontal line y = Vj, j = 0, 1, • • • ,N2, we obtain the new mesh 



such that 



'i+ij 



[^ij7 yij) 



Mj,(x, yij)dx = constant, i = 0,1, 



,iVi 



where M^ is the monitor function in the x coordinate direction. Note that 
only the first coordinates of the points have been changed (Fig. [1^). 

In the second stage the equi-distribution is performed in the vertical di- 
rection along the grid lines produced in the first stage. Since the grid lines 
are curved, the distribution is performed along the arc rather than the ver- 
tical coordinate. Denoting by s the arc-length variable for each vertical grid 



line i = 0, 1, • • • , A''i, the distance s 



«j' 



J 



0, 1,--- ,N2 from {x'-^,yio) to 
{x',-j,yij) along the vertical grid lines can be evaluated piecewise linearly by 

is^o = ) 

Sij = Sj(j-i) + \\{xij,yij) - (xj(j_i),yi(j_i))||. 





Figure 2: The subregion i^a in ^c is mapped into ftp. in il.p 



Having the values of the monitor function corresponding to Sij, i.e. the 
value of My at the points {x,-j,yij), the new e.d. mesh s^q,s^i, ■ ■ ■ ,s^j^^ is 
obtained by 






!{J + 1) 



My{x, y)ds = constant, i = 0, 1, ■ 



,iVi. 



The new values of s^ • can be used to generate the new positions of the points 
on the grid lines since the piecewise linear representation of the underlying 
arc is available (Fig. [lb). For more details see |16j . A similar procedure 
is performed in the third stage along the horizontal grid curves, with the 
monitor in the x coordinate direction again (Fig. [It). 
The mesh resulting from the above procedure forms quadrilaterals whose 
sides equi-distribute the grid lines in the two coordinate directions (see Fig. 



4 Adaptive nodes in a non-rectangular domain 

In this section we present a new technique for generating adaptive nodes in 
a simply connected domain, Qp, bounded by a closed curve T. For ease of 
explanation, we assume that T can be given by the parametric equations 

T:x = gi{e), y = g2{e), 0< ^ < 27r, 

where 9 is the angle used in the polar coordinate system and measured in 
the conventional anti-clockwise direction. The key to our proposed method 
is to make use of a rectangular transformation and to produce an adaptive 
mesh in the rectangle. We introduce a mapping 

V' : Oc ^ Op (3) 

where 

n^ = {{e,r),0 <9 <2Tr, < r < 1} 

is a rectangle in the cartesian coordinate system (9, r), referred to as compu- 
tational domain, Qp is the physical domain and the transformation is defined 
by the following functions 

x = rgi{9), ,. 

y = rg2{9). ^' 

One can easily see that the above transformation corresponds the line r = 1 
in ilc to the boundary T in Vtp. In addition, any rectangular subregion Qa 
in Oc is mapped into a subregion Qp. in Op as shown in Fig. [21 
The method described in section [3] and used to produce grid points in a rect- 
angle, is now applicable to the computational domain Vtc- Before explaining 
the adapting technique, we propose a method to construct a uniform mesh 
in Op by the use of the transformation in ([3]) . 

We start with a uniform mesh in Oc in the cartesian coordinate system 
9r as follows, 

6*0 < 6*1 < • • • < 9n-i = 27r - he, 

hr = ri <r2< ■■■ < tm-i <rm = l, 



where 



n 
and 



2tt 
he = — , 9j = jhe, j = 0, 1, . . . , n - 1, 
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Figure 3: The mesh points transformed from f^c into a circle before and 
after refinement are displayed in Figures (a) and (b) respectively. 

Using equations ^ the mesh points (6j,ri) are transformed to the points 
{xi,yj) in the xy coordinate system, as depicted in Fig. [3l Unlike the 
mesh in il.c, the new mesh in il.p is non-uniform, since each line r = r^ in 
Qc corresponds to a closed curve in fip with the same number of nodes, 
n. Therefore, the mesh will be clustered around the center, especially, for 
small values of r. This feature affects the quality of the grid points. To 
avoid this difficulty, below we suggest a refinement process to construct a 
roughly uniform mesh referred to as a uniform mesh in this paper. Although 
there might be some easier way to produce a set of uniform mesh points in 
an irregular region, the following method will be utilized in the adapting 
technique later. 



4.1 A uniform mesh in Qp 

In order to obtain a uniform mesh in Qp, we now refine the above mesh 
points obtained by the transformation. The refinement process is performed 
by modifying the number of nodes on each closed curve, r = ri based on 
a suitable criteria. For instance, we suggest this number of points to be 
selected in such away that the distance between the adjacent points on each 
closed curve be the same as that on the boundary, we propose the following 
process which is based on the perimeter of the closed curves. 
The perimeter of the ith curve can be approximated by the perimeter of a 
circle whose radius is evaluated by the average length of the position vectors 



of the boundary points, i.e. 

n-l 



Hi / ^ Hij ) i i ; ■ 



n 

j=0 



where 



Rij = ya^l + yfj, Xij = riQiiOj), Vij = rig2{0j). 

The approximate number of nodes for the ith curve can be therefore obtained 
by 

ne. = ^, (5) 

where pi = 27ri?j is the perimeter of the circle and As is the minimum 
distance between the adjacent boundary nodes. Another difficulty may arise 
due to the value of n^. In fact, if the values of ng and n^ are arbitrarily 
selected, the points may be clustered along the lines 9 = 9j (see Fig. H]). To 
treat this issue we propose the value of rij. to be selected based on ne using 
a criteria similar to the above mentioned process. For instance, for a given 
value of uq, rir can be determined such that 

P ng 

where p is the approximate perimeter of the boundary and R is the average 
lengthes of the position vectors. A clustered and a refined mesh obtained 
by the above process are displayed in Figures [3^ and[3]3, respectively. 
The suitably selected ng. and n^ will be also used in adapting mesh points 
later. 

4.2 Adaptive nodes in non-rectangular regions 

We now present a method to produce adaptive nodes in the physical domain 
Qp. The key to this new approach is to make use of the adaptive mesh 
technique described in section [3] and the transformation ([3]). More precisely, 
some adaptive points are first prepared in the computational domain f^c 
and then transformed into the physical domain fip. Fig. [5] illustrates how, 
for instance, the third stage of adapting mesh in f^c is related to 0,p, i.e. 
equi-distributing along the grid lines for the coordinate 9. 

Although, the mesh points produced by the above method is somehow 
adaptive, the concentration of the points around the center is an disadvan- 
tage as previously discussed. To overcome this difficulty, below we propose 
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Figure 4: The points produced by an arbitrarily selected value of n^ in (a) 
and a suitably selected value of n^ based on relation ([6]) in (b) are displayed. 




Figure 5: The grid curves in Qp resulted from the transformation of the grid 
lines in Qc in the third stage of adapting nodes. 
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a refining mechanism in performing the 3-stage adaptive algorithm. 
We start with a uniform mesh in J^c and perform the the first two stages of 
the adapting method. More precisely, the equi-distribution in the horizontal 
direction and the vertical grid lines are accomplished (Fig. [1^ and [1)3). But, 
the final stage is done in a different manner. We suggest a combination of 
the third stage of the adapting technique with the refinement process dis- 
cussed in section HTTl to avoid clustering the mesh points. As noted before, 
the density of the mesh around the center can be avoided by modifying the 
number of nodes along the closed curves. Recalling the refinement process to 
obtain a uniform mesh, a suitable number of points for the ith closed curve, 
ng. was computed in ([5]). We can either use the same number of points or 
a new value of ng- based on the new perimeter of the curves resulted from 
the first two stages. 

Having obtained the horizontal grid lines by the two-stage procedure and 
found the ng. for each curve, we can equi-distribute ng- points along the 
horizontal grid curves used in the third stage of adapting method in Fig. [1]:. 
Fig. [6] illustrates the 3-stage adapting nodes in the rectangular domain. It 
is observed that in the third stage (Fig. [6b) a refined number of nodes, ng. 
are distributed along the horizontal grid curves. 

5 Meshless methods 

In this section we first introduce the RBFs and then describe their ap- 
plication to the numerical solution of PDEs based on collocation meshless 
method. 

5.1 Radial basis functions 

RBFs are known as the natural extensions of splines to multi-variate inter- 
polation. Suppose the set of points 

{xien\i = i,2,--- ,N}, 

is given, where Q is a bounded domain in i?". The radial function cp : il — > 
R is used to construct the approximate function 

N 



fc=i 



which interpolates an unknown function / whose values at {xj}^^ are 
known. ||.|| represents the Euclidean norm. The unknown coefficients ak 

11 
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Figure 6: The illustration of the 3-stage adapting mesh in f^c with a refined 
number of nodes along each grid cure in the horizontal direction in (c) is 
displayed. 
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are determined such that the following A^ interpolation conditions are sat- 
isfied, 

N 

f{xi) = s{xi) = ^ak^{\\xi -Xfcll), i = I,--- ,A^. 
fc=i 
There is a large class of interpolating RBFs |14] that can be used in mesh- 
less methods. These include the linear 1 + r, the polynomial Pk{r), the 
thin plate spline (TPS) r^logr, the Gaussian exp(— r^//3^), and the multi- 
quadrics y /3^ -|- r^ (with /? a constant parameter). In this paper we employ 
a generalized TPS, i.e. r^logr which is a particular case of r^'^logr (k=2). 
These RBFs, augmented by some polynomials, are known as the natural ex- 
tension of the cubic splines to the case of 2D. In fact, they are obtained by 
minimizing a H"^ semi-norm over all interpolants for which the semi-norm 
exists. The theoretical discussion of these RBFs has been presented in [9]. 

5.2 Collocation meshless method 

We describe the collocation method for a general case of PDEs in the form 

Lu = F, (7) 

where L = [Li,--- ,Ln]'^ represents a vector of linear operations and 
F = [fi,- ■ ■ , /at]^ denotes a vector containing the right hand sides of the 
equations. For instance, Poisson's equation with a Dirichlet boundary con- 
dition 

Au = /, in il, 

u = g, on do,, 

is a very simple case of equation ([7]) where L = [A, I]'^, F = [f,g]'^ and the 
operators A and I act on the domain il. and the boundary dQ respectively. 
The collocation method is simply to express the unknown function u in 
terms of the RBFs as 

N 

u{x) = ^ak(p{\\x-Xk\\), (8) 

fc=i 

and determine the unknowns ak in such a way that ([8]) satisfies equation ([7]) 
for all interpolation points. Substituting ([8|) in equation ([7|) and imposing 
the N essential conditions of the collocation method lead to a linear system 
of equations whose coefficient matrix consists of A^ row blocks, the entries 
of which are of the form 

A'^j = Lf,ip{\\x - Xj\\)\x=a:^, i = !,■■■ ,Nf„ j = !,■■■ ,N, 

13 



where A'^ indicates the number of nodes associated with the operator L^. 
The above cohocation method is referred to as a non-symmetric coUocation 
inethod due to the non-symmetric coefficient matrices. The invertibihty of 
the coefficient matrix can not be guaranteed |10j, although in most cases a 
non-singular matrix is expected. To tackle this difficulty the symmetric col- 
location method, motivated by Hermitian interpolation, has been suggested. 
This method leads to symmetric matrices and proves the non-singularity of 
interpolation matrices at the expense of double acting the operators on the 
RBFs |llj . Since this work is not concerned with the singularity of the 
matrices the non-symmetric case will be implemented. Of course, the main 
technique discussed in this paper can be applied to the other case as well. 

6 Numerical Results 

We now examine the effectiveness of the mesh generation technique by ap- 
plying the collocation meshless method to some PDEs. In each case, the 
PDE is solved with some equally spaced nodes (roughly uniform) and adap- 
tive nodes generated by the new method and the results are compared. As 
previously noted, (/)(r) = r^logr, is used as a basis function. In each exam- 
ple, M test points, which do not coincide with the interpolation nodes, are 
randomly selected and a root mean square (RMS) error at these points is 
evaluated by 



RMS error 



M 



\ i=l 



where Uapr,i and Uex,i denote the approximate and exact values of u, respec- 
tively, at a test point i. 

Equations dH are used to evaluate the arc-length monitors 



Me = Jl + uj and M,. = y^l + u. 



respectively for equi-distributing in the 9 and r coordinate directions where 

du dx du dy du dx du dy 

dx dO dy d9 ' dx dr dy dr 

We consider the Poisson equation 

Au = f{x,y) in ri 
uix,y) = g{x,y), on dQ, 

14 



Table 1: Error values for Example [T] 



ng, Ur 


25,4 


31,5 


37,6 


43,7 


50, 8 


56,9 


63,10 


69,11 


n 


60 


86 


124 


162 


212 


266 


328 


398 



Adaptive 1.04E-3 8.47E-5 5.15E-5 2.31E-6 4.81E-6 5.14E-7 4.75E-7 7.67E-8 
Uniform 2.53E-3 9.24E-4 3.42E-4 1.36E-4 5.57E-5 2.37E-5 1.09E-5 4.80E-6 



Table 2: Error values for Example [2] 



ne, rir 


25,4 


31,5 


37,6 


43,7 


50, 8 


56,9 


63,10 


69,11 


n 


60 


86 


124 


162 


212 


266 


328 


398 



Adaptive 3.055E-4 3.76E-4 1.42-4 6.32E-5 1.91E-5 7.23E-6 2.50E-6 8.43E-7 
Uniform 2.89E-4 3.75E-4 2.26E-4 1.03E-4 4.32E-5 1.78E-5 7.83E-6 3.13E-6 



for different cases with various solutions and regions. 
Example 1 We solve equations ^ in the case of 

/(x, y) = 4e^'+2^' + 4(x2 + y2)e^'+J/' 

g{x,y) = e^ +y 

and 9il is an ellipse whose equation is x^ + 4^^ — 1 = 0. 
In this case the exact solution is given by u{x, y) = e^ ^^ . 

Example 2 

/(x, y) = -4e(4-^'-2^') + 4(x2 + y2)e(4-^'-s/') 

g{x,y) = e^^-^'-y") 
and dft is the same as that in Example {Ji 



The exact solution is given by u{x, y) = e 



15 



(4-^2 -y2) 



Table 3: Error values for Example [3] 



no, Ur 


25,4 


30,5 


45,8 


55,9 


65,11 


75,12 


n 


60 


78 


142 


262 


366 


460 



Adaptive 1.13E-3 3.80E-5 2.54E-5 3.21E-6 3.80E-7 1.81E-7 
Uniform 1.85E-3 7.62E-4 5.40E-5 2.27E-5 5.34E-6 2.46E-6 



Example 3 

g{x,y) = x^ + y2 
and the boundary is given by the parametric equations: 



X = cosOy^l — cos6/2 

y = sinQ^X — sinO /2 

The exact solution is given by u{x,y) = x^ + y^. 

The adaptive nodes generated for the above PDEs are displayed in Figures 
Uh, Ub and El respectively. One can expect, for the solution of Example (U 
the mesh points to be more dense close to the boundary and, for the solution 
of Example O the mesh points to be more dense around the center. These 
are exactly what we observe in the figures. Therefore, the adaptive nodes 
are in a good agreement with the solutions. A similar situation is observed 
in Fig. [8] for the solution of Example [31 

The numerical errors for the above three examples are listed, respec- 
tively, in Tables [H [2l and [3l where hq is the number of boundary points, n 
is the total number of collocation points and Ur is the appropriate number 
of divisions for the coordinate r based on relation (|6|) . 

In all the examples, the numerical results demonstrate a considerable reduc- 
tion in the error in the case of using the adapting nodes produced by the 
new adaptive method. In particular. Table [H shows that the error in the case 
of using 398 uniform nodes is nearly the same as that in the case of using 
212 adaptive nodes. Moreover, the error in the case of using 398 uniform 
nodes is 63 times larger than that in the case of using the same number of 
adaptive mesh points. 
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Figure 7: The adaptive nodes generated for the solutions of Examples [T] and 
[2] are displayed in (a) and (b) respectively. 
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Figure 8: The adaptive nodes generated for the solution of Example [3l 

7 Conclusion 

A new adaptive mesh technique was presented for irregular regions in two 
dimensional cases based on equi-distribution. The method was based on 
first, generating some adaptive nodes in a rectangle and then mapping the 
produced nodes to the physical domain using a suitable transformation. 

The new method was examined by considering some PDEs solved by a 
collocation meshless method and the results demonstrated considerable re- 
duction in the error values. Since the current work was not involved with 
constructing a monitor for the underlying method, a general monitor func- 
tion, arc-length, was employed to produce the adaptive nodes. Of course, 
using a suitable monitor for the underlying method could have resulted in 
more improvement in the numerical results 

The proposed method was implemented for 2D regions whose boundary 
was given by parametric equations in terms of polar coordinate system. 
For a more general case of irregular boundaries, we suggest a piecewise 
interpolation, in terms of some scattered points on the boundary, to make 
a set of parametric equations. Moreover, an extension of the new method 
to the case of 3D is currently under study and will be reported in a near 
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future. 
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